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COMPUTATION  OF  BESSEL  AND  NEUMANN  FUNCTIONS 

This  memorandum  describes  the  method  by  which  a program  was 
developed  for  evaluating  Bessel  and  Netimann  functions  of  Integral 
order  with  real  positive  arguments.  This  method  was  developed 
specifically  for  use  on  small  scale  computers  for  ^Ich  high  com- 
puting speed  can  be  obtained  If  excessive  precision  Is  avoided. 
The  program  allows  the  functions  to  be  computed  to  five  signifi- 
cant digits  using  six  digits  In  the  computation. 


The  work  described  In  this  memorandum  was  carried  out  to 
allow  computation  of  numerical  values  for  solutions  to  the  wave 
equation  in  cylindrical  coordinates  .j^This  problem  arose  in 
connection  with  a program,  of  work  In  Wogress  for  evalixatlng  the 
sound  pressure  field  near  a baffle- tram diicer  complex. 

Three  different  methods  are  used  t\evalilate  the  functions 
depending  on  the  value  of  the  order  and  thW  argument . The  range 
of  the  order  and  argument  for  each  of  the  thme  methods  is  shown 
in  Figure  1.  Note  that  the  three  regions  coiraletely  cover  all 
positive  Integral  values  of  order  and  all  positive  valties  of  the 
argument.  The  reasons  for  dividing  the  N-X  plane  into  the  regions 
shown  In  Figure  1 will  be  discussed  below. 

I.  BESSEL  FUNCTIONS 

A.  Region  1 

The  solution  of  Bessel's  differential  equation  of  the 
first  kind  is  usxially  referred  to  as  the  Bessel  function.  This 
function  may  be  expressed  as  a power  series, 
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(3) 

J^(x)  Is  the  Bessel  function  of  order  n and  argument  x.  The  dls.- 
cusslon  In  this  report  will  be  restricted  to  non-negative  Integral 
values  of  n and  real  non-negative  values  of  x.  The  power  series 
£>r  the  Bessel  function  Is  convergent  for  any  Integral  order  and 
any  argument.  The  usefulness  of  the  power  series  as  a computational 
algorithm  depends  on  the  behavior  of  the  terms  of  the  sum. 

Since  the  terms  of  the  sum  alternate  In  sign,  one  must  con- 
sider the  possibility  of  loss  of  significance  due  to  subtraction. 

To  examine  the  sijon  further,  consider  the  computation  of  J^(x), 


The  steps  In  evaluating  Jq(x)  for  x * 0,  1,  2,  6,  and  10 
are  shown  In  Table  1.  From  Table  1 It  can  be  seen  that  the  power 
series  converges  In  a few  terms  only  If  the  argument  Is  small. 

If  the  argument  Is  large  the  absolute  value  of  su^cesa^ve  terms 
diverges  until  the  term  number,  m,  exceeds  y ***  x * nl.  As 

can  be  seen  from  the  computation  of  J^(IO)  the  divergent  charac- 
teristics of  the  sum  leads  to  a subtractive  loss  of  significance. 
The  Insignificant  digits  of-Xhm  stans  are  circled.  Note  that  only 
three  significant  digits  are  obtained  In  the  computation  of  J^(IO) 
using  six  digits  of  significance  In  the  calculations.  Hence,  we 
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see  that  the  power  series  Is  a useful  computational  algorithm 
only  If  the  argument  Is  small  when  compared  to  the  order. 


To  obtain  a quantitative  expression  for  the  largest  value 
of  the  argument  \dilch  may  be  considered  "small,"  consider  equation 
(2)  for  m ■ 1, 


(~j) 

IHTTT  • 


If  the  ratio  of  to  Is  equated  to  -1,  the  resulting 

expression  can  be  solved  to  determine  the  value  of  the  argument 

for  which  T,  ■ -T^, 
i o' 

Xj^  - 2^n  + 1 . (6) 


If  the  argument  Is  less  than  or  equal  to  Xj^,  the  power 
aeries  may  be  used  to  compute  the  Bessel  function  since  the  ab- 
solute value  of  successive  terms  Is  non- Increasing.  If  k digits 
of  significance  are  used  In  computing  then  the  error  due 

to  subtraction  Is  approximately  • 10”^. 


B.  Region  2 

1.  General  Method 


On  Figure  1 regions  2 and  3 are  separated  by  the 
line  of  equal  order  and  argument.  For  values  of  order  greater 
than  two  the  Bessel  function  and  Its  derivative  can  be  evaluated 
for  equal  order  and  argument  uslng^. 


G.  N.  Watson,  A Treatise  on  the  Theory. of  Bessel  Functions. 
(Cambridge,  l^i!>o>,  /4o.  - — 
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If  J and  J'  are  known,  then  J ■,  and  J'  , can  be 

n n * n-1  n-1 

computed  by  the  recursion  relationships, 

(9) 


•>A-lW  - JnW  (10) 

In  theory  If  one  can  find  the  Bessel  function  and 
Its  derivative  for  any  order,  then  the  function  can  be  computed 
for  any  order  by  successive  application  of  Equations  (9)  and 
(10). 

In  practice  one  must  exercise  some  caution  In 
successive  application  of  recursion  formulae  as  subtractive  loss 
of  significance  can  occur.  The  loss  of  significance  Is  not 
severe  If  the  formulae  are  applied  In  region  2.  However,  se* 
vere  significance  loss  does  occur  If  these  formulae  are  applied 
with  Increasing  order  In  region  3. 

To  compute  a Bessel  function,  •l^(x),  in  region  2 
the  following  steps  are  carried  out: 

a.  The  argument  x Is  rounded  to  the  nearest 
integer,  x^. 

b.  (x^)  and  (x^)  are  computed  using  equations 
(7)  and  (8) 


Jn-l(*)  - JA(=‘)  + ; Jn(*). 
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c.  computed  by  successive 
application  of  the  recursion  formulae. 

d.  is  computed  using  a Taylor's  series  ex* 

pension  about  the  point  x^.  ^^^*0^  “^n^^o^ 

are  used  to  compute  the  coefficients  of  the 
2 

expansion. 

The  equations  for  the  Taylor's  srerles  expansion  are 
Included  here  for  completeness.  However,  the  derly/atlon  of  the 
equations  and  proof  of  convergence  are  net  Included  here. 


The  solutions  to  Bessel's  eqtiatlon,  S^(x),  may  be 
expressed  as  a Taylor's  series' expansion  about  the  point  x , 


Auatin  I,  Taxo* 


By  the  proper  choice  of  the  coefficients  A and  A, 

O L 

any  solution  of  Bessel’s  equation  can  be  generated.  Specifically 


1 . General  Method 

It  was  noted  earlier  that  successive  application  of 
the  recursion  fornnilae  in  region  3 with  increasing  order  leads  to 
severe  subtractive  loss  of  significance.  Hence,  the  basic  method 
which  was  applied  in  region  2 will  not  yield  satisfactory  results 
in  region  3. 

Earlier  it  was  shown  that  in  region  1 the  Bessel 
function  can  be  evaluated  using  the  power  series  expansion  about 
the  origin  shown  in  equation  (1).  Differentiation  of  equation  (1) 
yields 

“ /X\n+2m-l 

j:(x)  - T 
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To  compute  a Bessel  function,  Jjj(x),  in  region  3 the 
following  steps  are  carried  out: 

a.  Compute  ■ 2^n  + l'.  In  region  3,  x^  is  less 
than  X. 

b.  Compute  using  Equation  (1). 

c.  Compute  using  Eqxiation  (17). 

d.  Using  the  Taylor's  series  expansion  about  the 

point  Xq  compute  “^n^^o^  ^ere  x^  is 

either  x or  the  largest  value  of  the  argument 
which  may  be  used  without  severe  subtractive 
loss  of  significance  in  evaluating  the  Taylor's 
series  expansion. 

e.  If  it  is  not  possible  to  compute  Jj^(x)  in  step 
d,  the  values  computed  at  x^  are  used  to  repeat 
step  d until  it  is  possible  to  compute  JjjCx). 

Differentiation  of  Equation  (11)  yields  the  expression 
for  J*(x)  needed  in  step  d above. 


s;(x)  - 2 m • (18) 

m*l 


Equations  (11)  and  (18)  converge  in  a few  terms  if 
A is  small.  However,  if  A is  large.  Equations  (11)  and  (18)  dls 
play  slow  convergence  and  subtractive  loss  of  significance.  To 
avoid  this,  one  must  develop  some  method  of  choosing  an  upper 
limit  for  L.  Two  criteria  must  be  satisfied. 

1^1  1 0-9) 


and 
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(20) 


If  Equations  (19)  and  (20)  are  both  satisfied,  then 
Equations  (11)  and  (18)  will  converge  In  a few  terms  without 
severe  subtractive  loss  of  significance.* 


II.  NEUMANN  FUNCTIONS 
A.  Region  1 

The  solution  of  Bessel's  differential  equation  of  the 
second  kind  Is  usually  referred  to  as  the  Neumann  function.  This 
function  may  be  expressed  as  a series  expansion, 


n-1 


N„(x)  - 2[y  + tn(-^)J„(x)  - ^ Cn  - 


2m- n 


m"0 


n+2m 


2 k\l  +1) T + + 


m*0 


-1 1 

n + ml 


(21) 


where  N^(x)  “ The  Neumann  function  of  order  n and  argument  x; 

Y ■ Euler's  constant  ■ .5772157. 

Eqiiatlon  (21)  may  be  used  to  evaluate  the  Neumann 
function  In  region  1.  The  convergence  characteristics  and 


Equations  (19)  and  (20)  are  derived  from  the  recursion  Equation 
(15)  and  should  be  regarded  only  as  a sufficient  condition,  not 
a necessary  condition. 
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subtractive  loss  of  significance  are  very  similar  to  those  dis- 
cussed earlier  fbr.the  Bessel  fxinctlon. 

B.  Regions  2 and  3 

The  general  method  used  to  compute  Nexmoann  fvinctlons  In 
regions  2 and  3 Is  tlie  same  as  the  method  used  to  compute  Bessel 
functions  In  region  2. 

For  values  of  order  greater  than  two  the  Neumann  fvinctlon 
and  its  derivative  can  be  evaluated  for  equal  order  and  argument 
as. 


N„<n)  - 

.774759 

Fi  - ^ n 

1 .0101659 

fi 

(22) 

L 225n^. 

L 1462 Sn^J 

N;(n)  - 

,7ii6p  r. 

+ 23  1 

+ .154952  f. 

947  1 

(23) 

-;37r  [1 

3150n2j 

69300n2j  ‘ 

The  recursion  formulae  for  the  Neumann  functions  are 
Identical  to  those  for  the  Bessel  function,  so  equations  (9)  and 
(10)  may  be  restated  as, 

Nn-l(*>  " +5 

Cl<*)  - - N„(x).  (25) 

Successive  application  of  Equations  (24)  and  (25)  allows 
one  to  compute  N and  N'  for  any  Integral  order.  Severe  subtrac- 
tive loss  of  significance  does  not  occxjx  when  the  recursion 
formulae  ere  applied  to  the  Neumann  function  In  region  2 with  de- 
creasing order,  or  In  region  3 with  Increasing  order. 
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To  compute  a Neumann  function,  N^(x)  in  region  2 or  3 the 
following  steps  are  carried  out. 

1.  The  argument  is  rounded  to  the  nearest  integer,  x^, 

2.  N^  (x^)  and  N^  (x^)  are  computed  using  Eq\iations 
(22)  and '(23). 

3.  N^(Xq)  and  N^(Xq)  are  computed  by  successive  appli- 
cation of  the  recurrence  formulae 

4.  Nj}(x)  is  computed  using  the  Taylor's  series  expansion 
about  to  the  point  x^. 

III.  CONCLUSION 

The  method  developed  in  this  report  allows  one  to  compute 
the  Bessel  and  Neumann  functions  of  Integral  order  with  real 
positive  arguments.  Approximately  five  digits  of  accuracy  are 
obtained  using  six  digits  in  the  computation.  Subtractive  loss 
of  significance  has  been  minimized  in  the  method  developed  to 
avoid  time  consuming  multiple  precision  arithmetic  operations 
vdiich  would  otherwise  be  required  on  a small  scale  computer.  This 
allows  one  to  apply  a small  scale  computer  to  an  entire  class  of 
problems  which  were  previously  restricted  to  solution  on  a large 
scale  machine. 


